//===========================================================================
/*!
 *  \file Individual.cc

 *  \brief 
 *  
 *  \author  Nguyen Quoc Chinh
 *  \date    2007-01-25
  */
//===========================================================================
#include "Individual.h"

const char* Individual::Elements[]={"$",
	"H","He","Li","Be","B","C","N","O","F","Ne", //1st and 2nd round
	"Na","Mg","Al","Si","P","S","Cl","Ar", 	//3rd ro
	"K","Ca","Sc","Ti","V","Cr","Mn","Fe","Co","Ni","Cu","Zn","Ga","Ge","As","Se","Br","Kr"
	};

const char* Individual::AtomicNo[]={"0",
	"1","2","3","4","5","6","7","8","9","10",
	"11","12","13","14","15","16","17","18",
	"19","20","21","22","23","24","25","26","27","28","29","30","31","32","33","34","35","36"};

const double Individual::Mass[]={0,
	1.00783,4,6,8,10,12,14,15.99491,18,20,
	22.989,24.305,27,28,31,32,35.453,40,
	39,49,44,47.867,51,52,55,55.845,59,58.693,63.546,65.38,69.723,72.64,75,79,80,83.798};

Bond2Type_str Individual::bondOH={OXYGEN,HYDROGEN,0.2,1.2};
Bond2Type_str Individual::bondHO={HYDROGEN,OXYGEN,0.2,2.1};
Bond2Type_str Individual::bondOO={OXYGEN,OXYGEN,1.0,3.2};

Bond2Type_str Individual::bondFH={FLUORINE,HYDROGEN,0.2,1.4};
Bond2Type_str Individual::bondFF={FLUORINE,FLUORINE,1.0,3.2};

Bond2Type_str Individual::bondCO={CARBON,OXYGEN,1.0,1.5};
Bond2Type_str Individual::bondCH={CARBON,HYDROGEN,0.2,1.2};

Bond2Type_str Individual::HBondList[]={bondOH,bondFH,bondCH};
Bond2Type_str Individual::bondList[]={bondOH,bondOO,bondFH,bondFF,bondCO,bondCH};


Individual::Individual(){
	Individual::clear();
}

void Individual::clear(){	
	int i,j;
	nAtom=tag=0;
	energy=rmsGrad=0;
	for(i=0;i<DGR;i++) a[i]=alpha[i]=0;
	for(i=0;i<SIZE_MAX;i++){
		Nz[i]=0;
		for(j=0;j<DGR;j++){
			x[i*DGR+j]=force[i*DGR+j]=0;
		}
	}
	for(i=0;i<SIZE_MAX;i++) Nz[i]=0;
	for(i=0;i<TYPE_MAX;i++) nType[i]=0;
	for(i=0;i<12;i++) USRsig[i]=0;

	nFreq=0;
	iCorType=DIRECT;
}

Individual& Individual::operator=(const Individual& p) {
	int i;
	assert(p.nAtom < SIZE_MAX);
	for(i=0;i<DGR;i++){		a[i]=p.a[i];	alpha[i]=p.alpha[i];}
	nAtom = p.nAtom;	
	energy = p.energy;
	rmsGrad=p.rmsGrad;
	for(i=0;i<=nAtom*DGR; i++){	x[i]=p.x[i];	force[i]=p.force[i];	}	
	tag=p.tag;
	for(i=0;i<nAtom;i++) Nz[i]=p.Nz[i]; 
	nFreq=p.nFreq;
	for(i=0;i<nFreq;i++) freq[i]=p.freq[i];
	
	for(i=0;i<TYPE_MAX;i++) nType[i]=p.nType[i];
	
	for(i=0;i<12;i++) USRsig[i]=p.USRsig[i];
	iCorType=p.iCorType;
	//molList=p.molList;
    return *this;    // Return ref for multiple assignment
}//end operator=

void Individual::Migrate(Genes_str &p){
	int i;
	for(i=0;i<DGR;i++){
		p.a[i]=a[i];
		p.alpha[i]=alpha[i];
	}
	p.nAtom=nAtom ;	
	p.energy=energy;
	p.rmsGrad=rmsGrad;
	for(i=0;i <=nAtom*DGR; i++){	 p.x[i]=x[i];	 p.force[i]=force[i];}	
	p.tag=tag;
	for(i=0;i<p.nAtom;i++) p.Nz[i]=Nz[i];
	p.nFreq=nFreq;
	for(i=0;i<p.nFreq;i++) p.freq[i]=freq[i];
	for(i=0;i<TYPE_MAX;i++) p.nType[i]=nType[i];
	for(i=0;i<12;i++) p.USRsig[i]=USRsig[i];

	p.iCorType=iCorType;
}

void Individual::Get(Genes_str &p){
	int i;
	for(i=0;i<DGR;i++){
		a[i]=p.a[i];
		alpha[i]=p.alpha[i];
	}
	nAtom = p.nAtom;	
	energy = p.energy;
	rmsGrad=p.rmsGrad;	
	for(i=0;i <=nAtom*DGR; i++){
		x[i]=p.x[i];
		force[i]=p.force[i];
	}	
	tag=p.tag;
	for(i=0;i<nAtom;i++) Nz[i]=p.Nz[i];
	nFreq=p.nFreq;
	for(i=0;i<nFreq;i++) freq[i]=p.freq[i];
	
	for(i=0;i<TYPE_MAX;i++) nType[i]=p.nType[i];
	for(i=0;i<12;i++) USRsig[i]=p.USRsig[i];
	iCorType=p.iCorType;
}

void Individual:: AddStructure(Genes_str &p){	
	assert(nAtom+p.nAtom<SIZE_MAX);
	for(int i=0;i<p.nAtom;i++){
		Nz[nAtom+i]=p.Nz[i];
		nType[p.Nz[i]]++;
		for(int j=0;j<DGR;j++){
			x[(nAtom+i)*DGR+j]=p.x[i*DGR+j];			
		}		
	}
	nAtom+=p.nAtom;			
}

double Individual::torsion(int i1,int i2,int i3,int i4){
	double Q,P;
	
	double d12=distance(i1,i2);
	double d13=distance(i1,i3);
	double d14=distance(i1,i4);
	double d23=distance(i2,i3);
	double d24=distance(i2,i4);
	double d34=distance(i3,i4);

	P = SQR(d12) * ( SQR(d23)+SQR(d34)-SQR(d24)) +
		SQR(d23) * (-SQR(d23)+SQR(d34)+SQR(d24)) +
		SQR(d13) * ( SQR(d23)-SQR(d34)+SQR(d24)) -
		2 * SQR(d23) * SQR(d14);

	
	Q = (d12 + d23 + d13) * ( d12 + d23 - d13) *
	(d12 - d23 + d13) * (-d12 + d23 + d13 ) *
	(d23 + d34 + d24) * ( d23 + d34 - d24 ) *
	(d23 - d34 + d24) * (-d23 + d34 + d24 );
	
	//cout<<" P ="<<P<<" Q= "<<Q<<endl;
				
	return acos(P/sqrt(Q));
}

void Individual::RotateX(double theta){	
	double cost=cos(theta),sint=sin(theta);
	double xt,yt,zt;
	for(int i=0;i<nAtom;i++){
		xt=x[i*DGR]; yt=x[i*DGR+1]; zt=x[i*DGR+2];
		x[i*DGR+0]=  xt;		
		x[i*DGR+1]=  yt*cost + zt*sint;
		x[i*DGR+2]= -yt*sint + zt*cost;
	}
}

void Individual::RotateY(double theta){	
	double cost=cos(theta),sint=sin(theta);
	double xt,yt,zt;
	for(int i=0;i<nAtom;i++){
		xt=x[i*DGR]; yt=x[i*DGR+1]; zt=x[i*DGR+2];
		x[i*DGR+0]=  xt*cost +zt*sint;
		x[i*DGR+1]=  yt;
		x[i*DGR+2]= -xt*sint + zt*cost;
	}
}

void Individual::RotateZ(double theta){	
	double cost=cos(theta),sint=sin(theta);
	double xt,yt,zt;
	for(int i=0;i<nAtom;i++){
		xt=x[i*DGR]; yt=x[i*DGR+1]; zt=x[i*DGR+2];
		x[i*DGR+0]=  xt*cost + yt*sint;
		x[i*DGR+1]= -xt*sint + yt*cost;
		x[i*DGR+2]=  zt;
	}
}

void Individual::RotateAxis(double theta,double axis[]){	
	for(int i=0;i<nAtom;i++) RotateAxis_atom(i,theta,axis);
}

void Individual::Translate(double *vec,int direction){
	for(int i=0;i<nAtom;i++) Translate_atom(i,vec,direction);
}

void Individual::Transform(double D[][DGR]){
	for(int i=0;i<nAtom;i++) Transform_atom(i,D);
}

void Individual::Transform(double **D){
	for(int i=0;i<nAtom;i++) Transform_atom(i,D);
}

void Individual::RotateX_atom(int i,double theta){
	double cost=cos(theta),sint=sin(theta);
	double xt,yt,zt;
	
	xt=x[i*DGR]; yt=x[i*DGR+1]; zt=x[i*DGR+2];
	x[i*DGR+0]=  xt;
	x[i*DGR+1]=  yt*cost + zt*sint;
	x[i*DGR+2]= -yt*sint + zt*cost;
}

void Individual::RotateY_atom(int i,double theta){
	double cost=cos(theta),sint=sin(theta);
	double xt,yt,zt;
	
	xt=x[i*DGR]; yt=x[i*DGR+1]; zt=x[i*DGR+2];
	x[i*DGR+0]=  xt*cost +zt*sint;
	x[i*DGR+1]=  yt;
	x[i*DGR+2]= -xt*sint + zt*cost;
}

void Individual::RotateZ_atom(int i,double theta){
	double cost=cos(theta),sint=sin(theta);
	double xt,yt,zt;
	xt=x[i*DGR]; yt=x[i*DGR+1]; zt=x[i*DGR+2];
	x[i*DGR+0]=  xt*cost + yt*sint;
	x[i*DGR+1]= -xt*sint + yt*cost;
	x[i*DGR+2]=  zt;
}

void Individual::RotateAxis_atom(int i,double theta,double axis[]){
	double X,Y,Z,sum;
	sum=sqrt(SQR(axis[0])+SQR(axis[1])+SQR(axis[2]));
	X=axis[0]/sum;	Y=axis[1]/sum;	Z=axis[2]/sum;
	//construct rotation matrix
	double c,s,t;
	c=cos(theta);	s=sin(theta);	t=1-cos(theta);
	double R[4][4]={ t*X*X+c  , t*X*Y+s*Z, t*X*Z-s*Y,	0,
		t*X*Y-s*Z, t*Y*Y+c,   t*Y*Z+s*X,	0,
  t*X*Z+s*Y, t*Z*Y-s*X, t*Z*Z+c,		0,
  0,			0,			0,			1};
  int l,m,n;
  double xt[4],xtp[4];

  for(l=0;l<DGR;l++)	xt[l]=x[i*DGR+l];	xt[3]=1;
  MAT2_MUL_VEC(xtp,R,xt,4,4);
  for(l=0;l<DGR;l++){
	  x[i*DGR+l]=xtp[l];
  }
}

void Individual::RotateAxisWithPoint_atom(int i,double theta,double axis[],double point[]){
	Translate_atom(i,point,-1);
	RotateAxis_atom(i,theta,axis);
	Translate_atom(i,point,1);
}

void Individual::Translate_atom(int i,double *vec,int direction){
	x[i*DGR]+=direction*vec[0];
	x[i*DGR+1]+=direction*vec[1];
	x[i*DGR+2]+=direction*vec[2];
}

void Individual::Transform_atom(int i, double D[][DGR]){
	int l,m; double v[DGR];
	R(i,v);
	for(l=0;l<DGR;l++){
		x[i*DGR+l]=0;
		for(m=0;m<DGR;m++) x[i*DGR+l]+=D[l][m]*v[m];
	}
}

void Individual::Transform_atom(int i, double **D){
	int l,m; double v[DGR];
	R(i,v);
	for(l=0;l<DGR;l++){
		x[i*DGR+l]=0;
		for(m=0;m<DGR;m++) x[i*DGR+l]+=D[l][m]*v[m];
	}
}

void Individual::RotateAxis_mol(vector<int> &mk,double theta,double axis[]){
	double X,Y,Z,sum;
	sum=sqrt(SQR(axis[0])+SQR(axis[1])+SQR(axis[2]));
	X=axis[0]/sum;	Y=axis[1]/sum;	Z=axis[2]/sum;
	//construct rotation matrix
	double c,s,t;
	c=cos(theta);	s=sin(theta);	t=1-cos(theta);
	double R[4][4]={ t*X*X+c  , t*X*Y+s*Z, t*X*Z-s*Y,	0,
		t*X*Y-s*Z, t*Y*Y+c,   t*Y*Z+s*X,	0,
  t*X*Z+s*Y, t*Z*Y-s*X, t*Z*Z+c,		0,
  0,			0,			0,			1};
  int l,m,n;
  double xt[4],xtp[4];
  for(int i=0;i<mk.size();i++){
	  for(l=0;l<DGR;l++)	xt[l]=x[mk[i]*DGR+l];	xt[3]=1;
	  MAT2_MUL_VEC(xtp,R,xt,4,4);
	  for(l=0;l<DGR;l++)	x[mk[i]*DGR+l]=xtp[l];
  }
}

void Individual::Translate_mol(vector<int> &mk,double *vec,int direction){
	for(int i=0;i<mk.size();i++){
		x[mk[i]*DGR]+=direction*vec[0];
		x[mk[i]*DGR+1]+=direction*vec[1];
		x[mk[i]*DGR+2]+=direction*vec[2];
	}
}

void Individual::RotateAxisWithPoint_mol(vector<int> &mk,double theta,double axis[],double point[]){
	for(int i=0;i<mk.size();i++){
		RotateAxisWithPoint_atom(mk[i],theta,axis,point);
	}
}

void Individual::FixOrigin(){
	int i,k;
	//fix to origin
	double rc[DGR];
	CalcMassCenter(rc);
	for(i=0;i<nAtom;i++)
		for(k=0;k<DGR;k++) 
			x[i*DGR+k]-=rc[k];   
}

void Individual::StdOrient(){
	FixOrigin();
	double I[DGR][DGR],Ip[DGR]; 
	double X[DGR][DGR];
	int nrot;
	CalcInertiaTensor(I);
	//eigen33(I,Ip,X);
	jacobi33(I,Ip,X);
	//eigsrt(Ip,X,DGR);
	
	for(int i=0;i<DGR;i++,cout<<endl){
		cout<<Ip[i]<<" == ";
		for(int j=0;j<DGR;j++){	cout<<X[i][j]<<" "; }
	}
	
	Transform(X);
	jacobi33(I,Ip,X);
	for(int i=0;i<DGR;i++,cout<<endl){
		cout<<Ip[i]<<" == ";
		for(int j=0;j<DGR;j++){	cout<<X[i][j]<<" "; }
	}
	//FixOrigin();
}

double Individual::CalcUSRsimilar(Genes_str &p){
			double S =0;
			int i;
			for (i=0; i<12; i++) {
				S += fabs(USRsig[i] - p.USRsig[i]);
			}
			S /= 12.0;
			S = 1.0/(1.0+S);
			return S;
		}

double Individual::CalcRms(){
	rmsGrad=0;
	for(int i=0;i<nAtom*DGR;i++) rmsGrad+= SQR(force[i]);
	rmsGrad=sqrt(rmsGrad/(nAtom*DGR));
	return rmsGrad;
}

void Individual::CalcMassCenter(double rc[]){
	int i,k;
	for(k=0;k<DGR;k++) rc[k]=0;
	double mass=0;
	for(i=0;i<nAtom;i++){
		mass+=Mass[Nz[i]];
		for(k=0;k<DGR;k++)
			rc[k]+=Mass[Nz[i]]*x[i*DGR+k];
	}
	for(k=0;k<DGR;k++) rc[k]/=mass; 
}

void Individual::CalcInertiaTensor(double I[][DGR]){
	int l,m,i;
	MAT2_EQU_NUM(I,0,DGR,DGR);
	for(i=0;i<nAtom;i++){
		I[0][0]+= Mass[Nz[i]]*(SQR(x[i*DGR+1])+SQR(x[i*DGR+2]));
		I[1][1]+= Mass[Nz[i]]*(SQR(x[i*DGR+0])+SQR(x[i*DGR+2]));
		I[2][2]+= Mass[Nz[i]]*(SQR(x[i*DGR+0])+SQR(x[i*DGR+1]));

		I[0][1]+=-Mass[Nz[i]]*x[i*DGR+0]*x[i*DGR+1];
		I[0][2]+=-Mass[Nz[i]]*x[i*DGR+0]*x[i*DGR+2];
		I[1][2]+=-Mass[Nz[i]]*x[i*DGR+1]*x[i*DGR+2];
	}
	I[1][0]=I[0][1];
	I[2][1]=I[1][2];
	I[2][0]=I[0][2];
}

void Individual::CalcPrincipalMoments(double Ip[],double X[][DGR]){
	double I[DGR][DGR];
	int nrot;
	CalcInertiaTensor(I);
	//eigen33(I,DGR,Ip,X,&nrot);
	jacobi33(I,Ip,X);
}


bool Individual::CalcUSRsig(){
	//compute centroid (C)
	XYZ centroid;
	centroid.x = 0;
	centroid.y = 0;
	centroid.z = 0;
	int i;
		//clear USRsig
	for (i=0; i<12; i++){	USRsig[i] = 0;}
	
	if(nAtom==2){ USRsig[0]=distance(0,1); return true;}
	
	for (i=0; i<nAtom; i++) {
		centroid.x += x[i*DGR+0];
		centroid.y += x[i*DGR+1];
		centroid.z += x[i*DGR+2];
	}
	centroid.x /= (double) nAtom;
	centroid.y /= (double) nAtom;
	centroid.z /= (double) nAtom;

	//find distances to the centroid
	vector <double> distToCentroid;
	double maxDist = -1;
	double minDist = 1e14;
	int minItem = 0;
	int maxItem = 0;
	for (i=0; i<nAtom; i++) {
		double d = sqrt(pow(x[i*DGR+0] - centroid.x,2.0) + pow(x[i*DGR+1] - centroid.y,2.0) + pow(x[i*DGR+2] - centroid.z,2.0)  );
		if (d < minDist) {
			minDist = d;
			minItem = i;
		}
		if (d > maxDist) {
			maxDist = d;
			maxItem = i;
		}
		distToCentroid.push_back(d);
	}

	//find oxy closest to centroid (A)
	//A is minItem
	XYZ A;
	A.x = x[minItem*DGR+0];
	A.y = x[minItem*DGR+1];
	A.z = x[minItem*DGR+2];

	vector <double> distToA;
	for (i=0; i<nAtom; i++) {
		double d = sqrt(pow(x[i*DGR+0] - A.x,2.0) + pow(x[i*DGR+1] - A.y,2.0) + pow(x[i*DGR+2] - A.z,2.0)  );
		distToA.push_back(d);
	}

	//find oxy furthest from the centroid (B)
	//B is maxItem and update maxitem for C
	XYZ B;
	B.x = x[maxItem*DGR+0];
	B.y = x[maxItem*DGR+1];
	B.z = x[maxItem*DGR+2];

	vector <double> distToB;
	maxDist = -1;
	for (i=0; i<nAtom; i++) {
		double d = sqrt(pow(x[i*DGR+0] - B.x,2.0) + pow(x[i*DGR+1] - B.y,2.0) + pow(x[i*DGR+2] - B.z,2.0)  );
		distToB.push_back(d);
		if (d > maxDist) {
			maxDist = d;
			maxItem = i;
		}
	}

	//find oxy furthest from B (C)
	//C is in maxItem now
	XYZ C;
	C.x = x[maxItem*DGR+0];
	C.y = x[maxItem*DGR+1];
	C.z = x[maxItem*DGR+2];

	vector <double> distToC;
	for (i=0; i<nAtom; i++) {
		double d = sqrt(pow(x[i*DGR+0] - C.x,2.0) + pow(x[i*DGR+1] - C.y,2.0) + pow(x[i*DGR+2] - C.z,2.0)  );
		distToC.push_back(d);
	}

	//compute USRsigs from each distribution
	

	//compute means of the distributions
	for (i=0; i<nAtom; i++) {
		USRsig[0] += distToCentroid[i];
		USRsig[3] += distToA[i];
		USRsig[6] += distToB[i];
		USRsig[9] += distToC[i];
	}
	USRsig[0] /= (double) nAtom;
	USRsig[3] /= (double) nAtom;
	USRsig[6] /= (double) nAtom;
	USRsig[9] /= (double) nAtom;

	//compute variances of the distributions
	for (i=0; i<nAtom; i++) {
		USRsig[1] += pow(distToCentroid[i] - USRsig[0], 2.0);
		USRsig[4] += pow(distToA[i] - USRsig[3], 2.0);
		USRsig[7] += pow(distToB[i] - USRsig[6], 2.0);
		USRsig[10] += pow(distToC[i] - USRsig[9], 2.0);
	}
	USRsig[1] /= (double) (nAtom-1);
	USRsig[4] /= (double) (nAtom-1);
	USRsig[7] /= (double) (nAtom-1);
	USRsig[10] /= (double) (nAtom-1);

	//compute skewness of the distributions
	for (i=0; i<nAtom; i++) {
		USRsig[2] += pow(distToCentroid[i] - USRsig[0], 3.0);
		USRsig[5] += pow(distToA[i] - USRsig[3], 3.0);
		USRsig[8] += pow(distToB[i] - USRsig[6], 3.0);
		USRsig[11] += pow(distToC[i] - USRsig[9], 3.0);
	}
	USRsig[2] /= (double) (nAtom-1);
	USRsig[5] /= (double) (nAtom-1);
	USRsig[8] /= (double) (nAtom-1);
	USRsig[11] /= (double) (nAtom-1);
	USRsig[2] /= (double) pow(sqrt(USRsig[1]),3.0);
	USRsig[5] /= (double) pow(sqrt(USRsig[4]),3.0);
	USRsig[8] /= (double) pow(sqrt(USRsig[7]),3.0);
	USRsig[11] /= (double) pow(sqrt(USRsig[10]),3.0);

	return true;
};

void Individual::CorrectOrder(){
	int i,j,l,tempI;
	double tempD;
	for(i=0;i<nAtom-1;i++)
		for(j=nAtom-1;j>i;j--)
			if(Nz[j]>Nz[j-1]){
				swapAtom(j,j-1);		
			}
}

bool Individual::bCorrect(string &error){ error="";
	int d[SIZE_MAX][SIZE_MAX];	
	for(int i=0;i<nAtom;i++){
		for(int j=0;j<nAtom;j++){
			d[i][j]=0;
			for(int k=0;k<6;k++)	
				if(bBond2(i,j,bondList[k])){ 	
					d[i][j]=d[j][i]=1;
				break;
				}
		}
		d[i][i]=1;
	}
	/*
	cout<<endl<<"  ";
	for(int i=0;i<nAtom;i++) cout<<Elements[Nz[i]]<<" ";
	cout<<endl;
	for(int i=0;i<nAtom;i++){
		cout<<Elements[Nz[i]]<<" ";
		for(int j=0;j<nAtom;j++) cout<<d[i][j]<<" ";
		cout<<endl;
	}
	*/
	return bConnectivity(d,nAtom);
}

bool Individual::bBond2(int i1,int i2,Bond2Type_str  &bondType){
	bool answer=false;
	//if(Nz[i1]>=Nz[i2]) swap(i1,i2);
	if((Nz[i1]==bondType.t1)&&(Nz[i2]==bondType.t2)){
		double d=distance(i1,i2);
		answer=((d>=bondType.dmin)&&(d<=bondType.dmax));			
	}	
	return answer;
}

bool Individual::bConnectivity(int d[][SIZE_MAX],int n){
	int i,j,l;
	for(i=0;i<n;i++){
		for(j=0;j<n;j++)
			if((i!=j)&&(d[i][j]==1)){//if connected,jump to jth row and flip
			for(l=0;l<n;l++) if(d[i][l]==1) d[j][l]=1;
			}
	}
	for(i=0;i<n;i++)
		for(j=0;j<n;j++) if(d[i][j]==0) return false;

	return true;
}


int Individual::getAtomicType(string &s){
	int i,m1;
	s=trim(s);
	string t="";t=s[0];	t=t+s[1];	t=trim(t);
	s=s.substr(2);
	for(i=1;i<TYPE_MAX;i++){
		/*
		m1=s.find(Elements[i]);
		if(m1!=string::npos){
			//s.erase(0,m1+1);
			return i;
		}
		*/
		if((t.compare(Elements[i])==0)||(t.compare(AtomicNo[i])==0))  return i;
	}
	return 0;
}

void Individual::TIP_convert(int ow,int hw1,int hw2,Individual &water,int imodel){
	int i;
	//get o and h1,h2 coordinates
	double o[DGR],h1[DGR],h2[DGR];
	for(int l=0;l<DGR;l++){
		o[l]=x[ow*DGR+l];
		h1[l]=x[hw1*DGR+l];
		h2[l]=x[hw2*DGR+l];
	}
		
	//constructing template molecule
	double tip3p[9] ={0,0,0,	0.9572, 0.0000, 0.0000,	-0.2400, 0.9266, 0.0000};
	double tip4p[12]={0,0,0,	0.5859, 0.7570, 0.0000,	 0.5859,-0.7570, 0.0000,	0.1500, 0.0000,0.0000};
	double tip5p[12]={0,0,0,	0.4100,-0.5800,-0.6500,	-0.9300,-0.0800, 0.1700,	-0.06,-0.08,0.1};
	
	water.nAtom=0;
	water.Nz[0]=8;water.Nz[1]=1;water.Nz[2]=1;water.Nz[3]=1;water.Nz[4]=1;
	switch(imodel){
		case TIP3P:	water.nAtom=3;break;
		case TIP4P:	water.nAtom=4;break;
		case TIP5P:	water.nAtom=5;break;
		default: water.nAtom=3; 
			for(i=0;i<DGR;i++){
				water.x[0*DGR+i]=o[i];
				water.x[1*DGR+i]=h1[i];
				water.x[2*DGR+i]=h2[i];
			}
			return;
		break;
	}
		
	for(i=0;i<water.nAtom*3;i++){
		switch(imodel){
			case TIP3P:	water.x[i]=tip3p[i];break;
			case TIP4P:	water.x[i]=tip4p[i];break;
			case TIP5P:	water.x[i]=tip5p[i];break;		
		}		
	}
		
	//rotate template to the orientation of input
	double mid_inp[DGR]; //coordinates of the middle point of h1 and h2 of input
	double h_inp[DGR]; //coordinates of vectors(h1 to h2) of input
	double mid_temp[DGR]; //cooridinates of the middle point of h1 and h2 of template	
	double h_temp[DGR]; //coordinates of vectors(h1 to h2) of template	
	
	//calculate the orientation of molecules and templates (from oxygen to middle point of h1 and h2)
	double l_inp=0,l_temp=0;
	double axis[DGR];
		
	for(i=0;i<DGR;i++){		
		mid_temp[i]=(water.x[1*DGR+i]+water.x[2*DGR+i])/2-water.x[0*DGR+i];					
		h_temp[i]=water.x[1*DGR+i]-water.x[0*DGR+i];
		mid_inp[i]=(h1[i]+h2[i])/2-o[i];	
		h_inp[i]=h1[i]-o[i];		
	}
	//calculate the plane vector of inp and template
	double plane_inp[DGR],plane_temp[DGR];//unit vector
	CROSSPRODUCT(plane_temp,h_temp,mid_temp);
	NORMALISE(plane_temp);
	CROSSPRODUCT(plane_inp,h_inp,mid_inp);
	NORMALISE(plane_inp);
	CROSSPRODUCT(axis,plane_temp,plane_inp);
	double theta=acos(DOTPRODUCT(plane_temp,plane_inp));
	//and then rotate template around z axis by -phi
	water.RotateAxis(-theta,axis);
	
	for(i=0;i<DGR;i++){
		mid_inp[i]=(h1[i]+h2[i])/2-o[i];	
		mid_temp[i]=(water.x[1*DGR+i]+water.x[2*DGR+i])/2-water.x[0*DGR+i];			
	}		
	
	//calculate the angle between mid_inp and mid_temp
	l_temp=sqrt(DOTPRODUCT(mid_temp,mid_temp));
	l_inp=sqrt(DOTPRODUCT(mid_inp,mid_inp));
	double phi=acos(DOTPRODUCT(mid_temp,mid_inp)/(l_inp*l_temp));
	//calculate the axis for rotating
	CROSSPRODUCT(axis,mid_temp,mid_inp);
	//and then rotate template around z axis by -phi
	water.RotateAxis(-phi,axis);
	
	
	//translate template to oxygen's position
	water.Translate(o);	
}

void Individual::Write(ostream &os,int format){
		switch(format){
		case F_XYZ: WriteXYZ(os);	break;	//for MOVIE MOL format
		case F_CAR: WriteCAR(os);	break;	//for Acrylis format
		case F_GRO: WriteGro(os);	break; //for GROMACS format
		case F_PDB: WritePDB(os);	break; //for protein data bank format		
		case F_LATTICE: WriteXYZLattice(os);	break;	//for MOVIE MOL format	
		case F_GAUOUT: WriteGauOut(os); break;
		case F_INTERNAL: WriteInternal(os); break;
	}
}

void Individual::WriteXYZ(ostream &os){
	if(nAtom>0){
		os.setf(ios::fixed);
		//line1		
		os<<nAtom<<" "<<endl;
		//line2
		os<<format("%d %1.8lf %1.0e")%tag%energy%rmsGrad;
		os<<" %";//for(int i=0;i<SIG_MAX;i++) os<<" "<<sig[i];
		if(nFreq){
			os<<" @FREQ "<<nFreq;
			for(int i=0;i<nFreq;i++) os<<" "<<setprecision(4)<<freq[i];		
		}
		os<<" @IN"<<endl;
		//line3 and so on
		for(int i=0;i<nAtom;i++) {
			os<<Elements[Nz[i]]<<" ";
			//os<<Nz[i]<<" ";
			os<<setprecision(12)<<x[i*DGR+0]<<" "<<x[i*DGR+1]<<" "<<x[i*DGR+2];
			os<<setprecision(12)<<" "<<force[i*DGR+0]<<" "<<force[i*DGR+1]<<" "<<force[i*DGR+2];
			os<<endl;
		}
	}	
}

void Individual::WriteXYZLattice(ostream &os){
	if(nAtom>0){
		os.setf(ios::fixed);
		//line1		
		os<<nAtom<<" "<<a[0]<<" "<<a[1]<<" "<<a[2]<<" ";
		os<<alpha[0]*180.0/M_PI<<" "<<alpha[1]*180.0/M_PI<<" "<<alpha[2]*180.0/M_PI<<"\n";
		//line2
		os<<format("%d %1.8lf %1.0e")%tag%energy%rmsGrad;
		os<<" %";//for(int i=0;i<SIG_MAX;i++) os<<" "<<sig[i];
		if(nFreq){
			os<<" @FREQ "<<nFreq;
			for(int i=0;i<nFreq;i++) os<<" "<<setprecision(3)<<freq[i];		
		}
		os<<" @IN"<<endl;
		//line3 and so on
		for(int i=0;i<nAtom;i++) {
			os<<Elements[Nz[i]]<<" ";
			//os<<Nz[i]<<" ";
			os<<setprecision(12)<<x[i*DGR+0]<<" "<<x[i*DGR+1]<<" "<<x[i*DGR+2];
			//os<<setprecision(12)<<" "<<force[i*DGR+0]<<" "<<force[i*DGR+1]<<" "<<force[i*DGR+2];
			os<<endl;
		}
	}	
}

void Individual::WriteGro(ostream &os){
}

void Individual::WritePDB(ostream &os){
}

void Individual::WriteInternal(ostream &os){
	if(nAtom>=1){
		os.setf(ios::fixed);
		os<<nAtom<<endl;
		os<<format("%d %1.8lf %1.0e")%tag%energy%rmsGrad;
		os<<" %";//for(int i=0;i<SIG_MAX;i++) os<<" "<<sig[i];
		if(nFreq){
			os<<" @FREQ "<<nFreq;
			for(int i=0;i<nFreq;i++) os<<" "<<setprecision(4)<<freq[i];		
		}
		os<<" @IN"<<endl;
		//cout<<Elements[Nz[0]]<<endl;
	
		//if(nAtom>=2) os<<Elements[Nz[1]]<<" 1 "<<distance(1,0)<<endl;
		//if(nAtom>=3) os<<Elements[Nz[2]]<<" "<<distance(2,1)<<" "<<angle_dgr(2,1,0)<<endl;
		for(int i=0;i<nAtom;i++){
			double dist,distmin;
						
			int j1min=-1; distmin=1e9;
			for(int j=i-1;j>=0;j--){
				dist=distance(j,i);		if(dist<distmin){ distmin=dist; j1min=j;};
			}
			
			int j2min=-1;	distmin=1e9;
			for(int j=i-1;j>=0;j--) if(j!=j1min){
				dist=distance(j,j1min);		if(dist<distmin){ distmin=dist; j2min=j;};
			}
			
			int j3min=-1;	distmin=1e9;
			for(int j=i-1;j>=0;j--) if((j!=j1min)&&(j!=j2min)){
				dist=distance(j,j2min);		if(dist<distmin){ distmin=dist; j3min=j;};
			}
			
			os<<Elements[Nz[i]];
			if(j1min!=-1) os<<format(" %3d %8.3lf ")%(j1min+1)%distance(i,j1min);
			if(j2min!=-1) os<<format(" %3d %8.3lf ")%(j2min+1)%angle_dgr(j1min,i,j2min);
			if(j3min!=-1) os<<format(" %3d %8.3lf ")%(j3min+1)%torsion_dgr(i,j1min,j2min,j3min);
			os<<endl;
		}
	}
}

int Individual::Read(istream &is,int format){
	switch(format){
		case F_XYZ: ReadXYZ(is);	break;	//for xyz format	
		case F_LATTICE: ReadXYZ(is);	break;	//for xyzlattice format	
		case F_GRO: ReadGro(is);	break; //for GROMACS format
		case F_PDB: ReadPDB(is);	break; //for protein data bank format
		case F_CAR: ReadCAR(is);	break; //for Acrylis format
		case F_GAUOUT: ReadGauOut(is);	break; //for Gaussian output format
		case F_GAUOUT_C: ReadGauOut_C(is); break;
		case F_VASPOUT: ReadVASPOut(is); break; //for VASP output format
	}
	return nAtom;
}

int Individual::ReadXYZ(istream &is){
	clear();
	string s;
	getline(is,s);
	//sscanf(s.c_str(),"%d",&nAtom);	
	sscanf(s.c_str(),"%d %lf %lf %lf %lf %lf %lf ",&nAtom,&a[0],&a[1],&a[2],&alpha[0],&alpha[1],&alpha[2]);
	
	for(int i=0;i<DGR;i++) alpha[i]*=M_PI/180.0;
	
	//cout<<s<<"\n nsdsf="<<nAtom<<endl;
	getline(is,s);
	sscanf(s.c_str(),"%lf",&energy);	
	
	if(s.find("@IN")!=string::npos){ //for my format
		sscanf(s.c_str(),"%d %lf %lf",&tag,&energy,&rmsGrad);		
		if(s.find("@FREQ")!=string::npos){ // if it has frequencies
				sscanf(s.c_str(),"%d %lf %lf @FREQ",&tag,&energy,&rmsGrad);				
				s.erase(0,s.find("@FREQ")+6);
				stringstream ss; ss<<s<<" ";
				
				ss>>nFreq;				
				//cout<<s<<" == "<<nFreq<<endl;
				for(int i=0;i<nFreq;i++) ss>>freq[i];			
		}
	}else if(s.find("TTM21")!=string::npos){ //for Harold format
			double eTTM2F,gradTTM2F,eOSS2,gradOSS2;
			sscanf(s.c_str(),"TTM21: %lf RMS: %lf  OSS2 : %lf RMS: %lf",&eTTM2F,&gradTTM2F,&eOSS2,&gradOSS2);
			energy=eOSS2;
			rmsGrad=gradOSS2;
		}
	
	int i,k;
	
	for(i=0;i<nAtom;i++){
		getline(is,s);		s+=" 0 0 0 ";
		Nz[i]=k=getAtomicType(s);
		nType[k]++;
		if(k>0){			
			sscanf(s.c_str(),"%lf %lf %lf %lf %lf %lf",&x[i*DGR+0],&x[i*DGR+1],&x[i*DGR+2],&force[i*DGR+0],&force[i*DGR+1],&force[i*DGR+2]);		
		}
		//else cout<<" Strange element ! ";
	}

	return 0;
}

int Individual::ReadVASPOut(istream &is){
	clear();
	
	return 0;
}

int Individual::ReadGro(istream &is){
	clear();	
	return 0;
}

int Individual::ReadPDB(istream &is){
	clear();	
	string s,temp;	
	int i=0,m,k;	
	while(!is.eof()){		
		getline(is,s);
		if((s.find("CRYST1"))!=string::npos){
			sscanf(s.c_str(),"CRYST1  %lf %lf %lf %lf %lf %lf ",&a[0],&a[1],&a[2],&alpha[0],&alpha[1],&alpha[2]);
			for(m=0;m<DGR;m++) alpha[i]*=M_PI/180.0;
		}

		if((s.find("ATOM"))!=string::npos){
			temp=s;
			k=20;temp.erase(0,k);			
			sscanf(temp.c_str(),"%d       %lf   %lf   %lf",&k,&x[i*DGR+0],&x[i*DGR+1],&x[i*DGR+2]);
			//for(m=0;m<DGR;m++) x[i*DGR+m]/=10.0;
			m=0; temp=s; temp.erase(0,12); temp.erase(3);temp=trim(temp);			
			k=getAtomicType(temp);
			//cout<<temp<<"  & "<<nAtom<<" & "<<k<<endl;
			if(k){		nAtom++;	nType[k]++;	Nz[i]=k;	i++;	}
		}
	}
	return 0;
}

int Individual::ReadGauOut(istream &is){
	clear();
	int flag=0,i,t1,pos1,pos2;
	string s,tmp;
	while(!is.eof()){
		getline(is,s); tmp=trim(s);
		
		//read atomic number and coordinates
		
		if(((s.find("Input orientation",0)!= string::npos)||(s.find("Standard orientation",0)!= string::npos))
				   &&(flag==0)){
			getline(is,s);getline(is,s);getline(is,s);getline(is,s);
			nAtom=i=0;
			//cout<<s<<endl;
			while(getline(is,s)){
				if(s[1]=='-') break;
				sscanf(s.c_str(), "%d %d %d %lf %lf %lf", &t1,&Nz[i],&t1,&x[i*DGR],&x[i*DGR+1],&x[i*DGR+2]);
				i++;
			}
			nAtom=i;
			
			flag|=1;//indicate coordinates have been collected
		}
			
		//read the energy
		if((s.find("E(RB+HF-LYP)",0)!=string::npos)&&(flag==1)){
			pos1=s.find("E(")+15;
			pos2=s.find("A.U.");
			sscanf(s.substr(pos1,pos2-pos1).c_str(),"%lf",&energy);
			flag|=2; //indicate energy has been collected
		}

		//read the RMS
		if((s.find("RMS     Force",0)!=string::npos)){
			pos1=25;			pos2=37;
			sscanf(s.substr(pos1,pos2-pos1).c_str(),"%lf",&rmsGrad);
		}

		
		if ((tmp.find("EUMP2",0) != string::npos)&&(flag==1)){
                // wow, energy
			tmp.replace(0,34,"");
			int dd=tmp.find("D",0);
			tmp.replace(dd,1,"E");
			sscanf(tmp.c_str(), "%lf", &energy);
			flag|=2; //indicate energy has been collected
		}

		if(flag==3){ 	flag=0;	}
		
		if(s.find("Frequencies --")!= string::npos){
			tmp=trim(s);
			//cout<<tmp<<endl;
			sscanf(tmp.c_str(), "Frequencies -- %lf %lf %lf",&(freq[nFreq]),&(freq[nFreq+1]),&(freq[nFreq+2]));				
			nFreq+=3;
		}		
		
		//if found the error in termination
		//if(s.find("Error termination")!=string::npos){ return -2;}
	}	
	return 0;
}


int Individual::ReadGauOut_C(istream &is){
	clear();
	int flag=0,i,t1,pos1,pos2;
	string s,tmp;
	while(!is.eof()){
		getline(is,s); tmp=trim(s);
		
		//read atomic number and coordinates
		
		if(((s.find("Input orientation",0)!= string::npos)||(s.find("Standard orientation",0)!= string::npos))
				   &&(flag==0)){
			getline(is,s);getline(is,s);getline(is,s);getline(is,s);
			nAtom=i=0;
			//cout<<s<<endl;
			while(getline(is,s)){
				if(s[1]=='-') break;
				sscanf(s.c_str(), "%d %d %d %lf %lf %lf", &t1,&Nz[i],&t1,&x[i*DGR],&x[i*DGR+1],&x[i*DGR+2]);
				i++;
			}
			nAtom=i;
			
			flag|=1;//indicate coordinates have been collected
		}
			
		//read the energy
		if((s.find("E(RB+HF-LYP)",0)!=string::npos)&&(flag==1)){
			pos1=s.find("E(")+15;
			pos2=s.find("A.U.");
			sscanf(s.substr(pos1,pos2-pos1).c_str(),"%lf",&energy);
			flag|=2; //indicate energy has been collected
		}

		//read the RMS
		if((s.find("RMS     Force",0)!=string::npos)){
			pos1=25;			pos2=37;
			sscanf(s.substr(pos1,pos2-pos1).c_str(),"%lf",&rmsGrad);
		}

		
		if ((tmp.find("EUMP2",0) != string::npos)&&(flag==1)){
                // wow, energy
			tmp.replace(0,34,"");
			int dd=tmp.find("D",0);
			tmp.replace(dd,1,"E");
			sscanf(tmp.c_str(), "%lf", &energy);
			flag|=2; //indicate energy has been collected
		}

		if(flag==3){ 	flag=0;	}
		
		if(s.find("Frequencies --")!= string::npos){
			tmp=trim(s);
			//cout<<tmp<<endl;
			if(nAtom==2){
				sscanf(tmp.c_str(), "Frequencies -- %lf",&(freq[nFreq]));				
				nFreq++;
			}else{
				sscanf(tmp.c_str(), "Frequencies -- %lf %lf %lf",&(freq[nFreq]),&(freq[nFreq+1]),&(freq[nFreq+2]));				
				nFreq+=3;
			}			
		}		
		
		//if found the error in termination
		if((s.find("Error termination")!=string::npos)||(s.find("cpu")!=string::npos)){ return -2;}
	}	
	return 0;
}


//============ Supporting functions ============/

string Individual::trim(string s){
	while (s.length() > 0 && s[0] == ' ') s = s.substr(1);
	while (s.length() > 0 && s[0] == '\t') s = s.substr(1);
	while (s.length() > 0 && s[s.length()-1] == ' ') s = s.substr(0, s.length() - 1);
	return s;
}

#define ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);\
	a[k][l]=h+s*(g-h*tau);
int Individual::jacobi33(double a[][3],double d[], double v[][3])
//Computes all eigenvalues and eigenvectors of a real symmetric matrix a[0..n-1][0..n-1]. On
//output, elements of a above the diagonal are destroyed. d[0..n-1] returns the eigenvalues of a.
//v[0..n-1][0..n-1] is a matrix whose columns contain, on output, the normalized eigenvectors of
//a. nrot returns the number of Jacobi rotations that were required.
{
	int j,iq,ip,i;
	double tresh,theta,tau,t,sm,s,h,g,c,*b,*z;
	int n=3;

	b=new double[n];
	z=new double[n];
	for (ip=0;ip<n;ip++) {						//Initialize to the identity matrix.
		for (iq=0;iq<n;iq++) v[ip][iq]=0.0;
		v[ip][ip]=1.0;
	}
	for (ip=0;ip<n;ip++) {						//Initialize b and d to the diagonal
												//of a. 
		b[ip]=d[ip]=a[ip][ip];
		z[ip]=0.0;									//This vector will accumulate terms
													//of the form tapq as in equa-
													//tion (11.1.14).
	}
	int nrot=0;
	for (i=1;i<=100;i++) {
		sm=0.0;
		for (ip=0;ip<n-1;ip++) {						//Sum o-diagonal elements.
			for (iq=ip+1;iq<n;iq++)
				sm += fabs(a[ip][iq]);
		}
		if (sm == 0.0) { //The normal return, which relies on quadratic convergence to machine under ow.
			//eigsrt(d,v,n);
			delete[] b;
			delete[] z;
			return nrot;
		}
		if (i < 4)
			tresh=0.2*sm/(n*n);								// ...on the rst three sweeps.
		else
			tresh=0.0;										//...thereafter.
		for (ip=0;ip<n-1;ip++) {
			for (iq=ip+1;iq<n;iq++) {
				g=100.0*fabs(a[ip][iq]);
				//After four sweeps, skip the rotation if the o-diagonal element is small.
				if (i > 4 && (double)(fabs(d[ip])+g) == (double)fabs(d[ip])
								&& (double)(fabs(d[iq])+g) == (double)fabs(d[iq]))
					a[ip][iq]=0.0;
				else if (fabs(a[ip][iq]) > tresh) {
					h=d[iq]-d[ip];
					if ((double)(fabs(h)+g) == (double)fabs(h))
						t=(a[ip][iq])/h;				//t = 1=(2.)
					else {
						theta=0.5*h/(a[ip][iq]);		//Equation (11.1.10).
						t=1.0/(fabs(theta)+sqrt(1.0+theta*theta));
						if (theta < 0.0) t = -t;
					}
					c=1.0/sqrt(1+t*t);
					s=t*c;
					tau=s/(1.0+c);
					h=t*a[ip][iq];
					z[ip] -= h;
					z[iq] += h;
					d[ip] -= h;
					d[iq] += h;
					a[ip][iq]=0.0;
					for (j=0;j<ip-1;j++) {					//Case of rotations 1 . j < p.
						ROTATE(a,j,ip,j,iq)
					}
					for (j=ip+1;j<iq-1;j++) {					//Case of rotations p < j < q.
						ROTATE(a,ip,j,j,iq)
					}
					for (j=iq+1;j<n;j++) {						//Case of rotations q < j . n.
						ROTATE(a,ip,j,iq,j)
					}
					for (j=0;j<n;j++) {
						ROTATE(v,j,ip,j,iq)
					}
					++nrot;
				}
			}
		}
		for (ip=0;ip<n;ip++) {
			b[ip] += z[ip];
			d[ip]=b[ip];								//Update d with the sum of tapq,
			z[ip]=0.0;									//and reinitialize z.
		}
	}
	return nrot;
	//nrerror("Too many iterations in routine jacobi");
}


#undef ROTATE